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A model is tested to rapidly evaluate the vibrational properties of alloys with site disorder. It 
is shown that length-dependent transferable force constants exist, and can be used to accurately 
predict the vibrational entropy of substitutionally ordered and disordered structures in Au-Cu, Au- 
Pd, and Cu-Pd. For each relevant force constant, a length- dependent function is determined and 
fitted to force constants obtained from first-principles pseudopotential calculations. We show that 
these transferable force constants can accurately predict vibrational entropies of Ll2-ordered and 
disordered phases in CusAu, AuaPd, PdsAu, CusPd, and PdaAu. In addition, we calculate the 
vibrational entropy difference between Ll2-ordered and disordered phases of AuaCu and CusPt. 

I. INTRODUCTION 

In the last decade, a clear prescription has emerged to obtain phase diagrams of materials from first-principles. When 
two compounds mix, some amount of site disorder occurs. First-principles alloy theory has been mostly preoccupied 
with defining appropriate energy models for such partially disordered systems, and obtaining the parameters for them 
from increasingly accurate Density Function Theory (DFT) methods. In particular, the cluster expansion approach ^ 
has been highly successful, as it allows one to parameterize the energy of systems with disorder on a fixed framework 
of sites. Many binary phase diagrams have been calculated in this way . 

Most of these first-principles phase diagram calculations do not include vibrational entropy effects. The assumption 
that vibrational entropy differences between phases are small, and can therefore be neglected, has been recently 
examined experimentally ^^^^ and theoretically 12-18 jj^ AI2CU vibrational entropy contributions were shown 
to be essential in reversing the stability of and 0' phases at finite temperatures. In the Al-Sc system, vibrational 
entropy contributions were shown to increase the solubility limits 27-fold i^'^'^. 

For alloys with partial site disorder, vibrational entropy in first-principles phase diagrams is computationally inten- 
sive to determine. In conventional first-principles alloy theory, a cluster expansion has to be fit to the ground state 
energies of a large number of A-B ordered states (for a binary alloy). Vibrational effects can be formally included 
by fitting to vibrational free energies, rather than ground state energies This requires that one determines the 
phonon spectrum for many configurations. Unfortunately, this is computationally expensive. 

One idea that has been proposed to alleviate this problem is the use of force constants that can be transferred 
between different environments. The existence of transferable force constants would greatly reduce the number 
of calculations required to determine vibrational entropies from first-principles: currently, almost all of the work in 
calculating vibrational entropies is involved with the numerical determination of force constants for every environment 
in a large number of structures. 

The idea of using transferable force constants to reduce computational expense is tremendously appealing. Despite 
this, until recently, there was little clear evidence supporting the validity of their use. For instance, experimental 
data in FcaAl shows that force constants can have a clear configurational dependence Strictly speaking, it is 
not possible to rigorously define configuration-independent force constant matrices as the form of the force constant 
matrix depends on the symmetry of the configuration Nevertheless, configuration-independent transferable force 
constants have been used in computations to study segregation in NisAl and elastic properties in Li-Al . In Li-Al, 
these force constants predicted elastic constants well, but reproduced vibrational entropies poorly. Transferable force 
constants have also been explored in oxides and semiconductors In these materials, transferability is examined 
in the short-ranged force constants, which are obtained after subtracting out an analytical long-ranged dipole-dipole 
contribution. None of these studies examined the use of transferable force constants to calculate vibrational entropy. 
This is a challenging issue, as vibrational entropy differences between structures are a fraction of the total vibrational 
entropy and hence, high relative accuracy is required. 

Previous attempts at defining transferable force constant matrices defined force constant matrices for a given 
pair type. Recently, it was shown that high-accuracy in predicting the vibrational entropy could be achieved by using 
bond length-dependent transferable force constants 12,26 rpj^-^ gjjj^pjg bond length-dependent force constant model 
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was shown to work well in predicting vibrational entropies in NisAl and PdsV 1^,13 q^i^^ force constants were found to 
be transferable between different FCC structures at a single composition; their transferability between configurations 
with different atomic ratios or different chemical systems was untested. 

In this paper, we evaluate the accuracy of the length-dependent transferable force constants approach in calculating 
vibrational entropies in the Au-Cu, Au-Pd, and Cu-Pd systems, and also investigate the transferability of such relations 
between different chemical systems. Our results show that in the Au-Cu, Au-Pd, and Cu-Pd systems, transferable 
force constants can be defined that depend on length and pair type. We show that these force constants can be 
exchanged between different structures and chemical systems, while still retaining good accuracy for calculating the 
vibrational entropy of intermetallic structures. 



II. METHODOLOGY 



A supercell method was used to calculate force constants 12,26-29 r^j^-g jjjq^JjqjJ consists of creating a series of 

supcrcells with well-chosen displacements, calculating the resulting forces on atoms, then using a least-squares fit 
to obtain the force constants. All total energy and force calculations were performed using the Vienna ab initio 
Simulation Package (VASP) 30,3i_ -^yi^f-i^ implements Blochl's projector augmented wave (PAW) approach '^^■^^ within 
the local density approximation (LDA). The calculations used energy cutoffs of 300-400 eV and ^150-200 unique 

fc-points per 4 atoms. This resulted in energies converged to ~5 meV, and forces converged to 1%. Special care was 
taken to make sure forces on the initial structure vanished before supercell configurations were generated for the force 

constant calculation. 

The superstructures and their perturbations were determined using an efficient algorithm which generates 
supercells with perturbed atoms and their corresponding set of linear equations. Third-order anharmonic terms, 

which are present when positive and negative displacements are not symmetrically equivalent, were eliminated by 
subtracting forces taken from calculations of perturbations of equal magnitude but opposite sign Vibrational 

entropy calculations were converged to within 0.015 ks- this value is obtained by calculating the configuration- 

00 

dependent part of the high-temperature limit of the vibrational entropy per atom, Syib = —ks J g{co)lnco dco, for 



increasing force constant range until S'l.ibchanged by less than 0.015 ks- Typically, this requires ~4-5 neighbors for 
LI2 structures and 2 neighbors for SQS-8 structures. 

The above procedure was applied to calculate vibrational entropies of Ll2-ordered and disordered phases in the 
AusCu, CusAu, AuaPd, PdsAu, CuaPd, PdsAu, and CusPt systems. In addition, we calculated vibrational entropies 
of Ll2-ordered Ag3Au and AusAg. The disordered structure was approximated by an 8 atom special quasirandom 
structure (SQS) These structures have been shown to give the best possible approximation to a disordered 
structure, within a given number of sites. In NiaAl, Morgan et. al. tested the convergence of vibrational properties 
with SQS size, and found that a SQS with 8 atoms is a good approximation of the disordered state The SQS-8 
used in this study has been used successfully in first-principles vibrational entropy calculations of disordered NisAl 
and PdsV ^^.la. 

To define transferable force constant matrices, some approximations are necessary. The form of a force constant 
matrix, and hence the number of non-zero force constants, depends on the symmetry of structure. In this work, we 
use a stretching-bending force constant model. In this model, the coordinate system of each force constant matrix 
is transformed, so that the 2-axis is aligned along the segment joining the two atoms in question. Two further 
approximations are necessary to obtain transferable force constants. First, bending terms arc averaged, so that they 
are orientation independent. Second, off-diagonal terms are constrained to be zero. The resulting force constant 
matrix has only two independent terms - a stretching term s, and a bending term b - and the form of the matrix can 
be written 



fb 0\ 

6 0. (1) 
\0 s J 

The stretching term (s) can be used as a qualitative measure of the bond strength. However, it is not possible to 
keep only the stretching terms (s) in the force constant matrix: doing so can result in large errors (0.2 ks) in the 
vibrational entropy 12,26 examined two types of errors associated with using stretching-bending force constant 
matrices: errors introduced in the force constants and errors introduced in the vibrational entropy. For all force 
constant matrices examined, the rms error associated with constraining the bending terms to be equal was 0.052 
eV/A? per bending term. The rms error associated with setting off-diagonal terms to zero was 0.009 eV/A^ per 
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ofF-diagonal term. Stretching terms contain no errors associated with the stretching-bending force constant model: 
they are directly obtained from the coordinate transformation, with no further approximations. 

The errors introduced into the vibrational entropy by using the stretching-bending force constant model are shown in 
Table 1. For 14 out of the 16 structures tested, Syib calculated using full-force constant matrices and Syib calculated 
using the stretching-bending force constant model agreed to within 0.01 ks- For the LI2 AuaCu structure, the 
stretching-bending force constant approximation failed when using nearest-neighbor force constants; for the SQS-8 
AusCu structure, the stretching-bending force constant approximation failed when using both nearest-neighbor and 
all-neighbor force constants (where "all-neighbor" is defined as the force constant range at which vibrational entropies 
were converged to within 0.015 ks)- In these structures, using stretching-bending force constant matrices resulted 
in a dynamically unstable structure. Tests of each approximation revealed that the negative phonon modes in the 
LI2 AuaCu structure arise because the off-diagonal terms to set to zero; negative phonon modes in the SQS-8 AuaCu 
arise due to both the neglect of the off-diagonal terms and the constraint on the bending terms. In previous work, the 
errors in vibrational entropy by using the simplified stretching-bending force constant model have been tested for five 
structures in the NiaAl and PdsV system 12, 13 those systems, the errors on vibrational entropy for all structures 
were less than 0.01 fcs/atom. 



III. RESULTS AND DISCUSSION 



Calculated vibrational entropies for all structures examined in the Ag-Ag, Au-Cu, Au-Pd, Cu-Pt, and Cu-Pd 
system are listed in Table 1. For the Cu-Au system, our results can be compared with previous experimental 
and theoretical work. We calculated Svib (Cu) = 4.87 fcs and Svib (^w) =2.85 ks, which yields for the formation 
entropiesS^°f,™(^l2 CusAu) = 0.10 ks and S'^°'''"(Ll2 AmCu) = 0.11 ks- Our calculated S'^°'''"(Ll2 Cu3Au)agTees 
well with previous experimental work, which found Sl°™{Ll2 Cu^Av) = 0.07 ± 0.03 ks ■ Our calculated 
formation entropies also agree well with previous theoretical work, which found S^"™ {LI2 Cu^Au) = 0.10 ks and 
Sia"'{Ll2 AusCu) = 0.14 ks The entropy change upon disordering in CusAu l^l:^s°rder^disorder (Cu^Au)) is 
calculated to be 0.07 ± 0.045 ks- This value is lower than the experimental value of 0.14 ± 0.05 ks ■ However, 
this calculation is in good agreement with previous theoretical work, in which /^^gorder^disorder ((7^3^^) = 0.O8 ks 
was obtained For AusCu, we calculate AS'°[^'^'=''^*«°'''^«'' (AusCu) = -0.01 ± 0.03 fcs, which is lower than the 
previously calculated AS'°[^'''='-^*"°'^'''='- (AusCu) = 0.05 fcs 

Calculated full force constant matrices were transformed to the stretching-bending force constant model, with the 
results plotted in Figures 1-6. The different bond lengths correspond to equilibrium bond lengths in various structures. 
For all pair types, the force constant stiffness for both stretching and bending terms decreases with increasing bond 
length. In particular, the stiffness of first nearest-neighbor force-constant matrices shows a strong dependence on 
bond length. This dependence of bond stiffness on bond length shows that one can not rely too much upon simple 
bond-counting arguments when predicting or explaining vibrational entropy differences between phases because 
the bond counting effect does not take into account bond-length changes upon disordering. Figures 1-6 show that 
small changes in bond length upon disordering can have a large effect on force constants and vibrational entropy. The 
relationship between atomic relaxations and vibrational entropy has been previously noted 12,16-18 

For a given bond type, stretching and bending force constants for all systems and structures are found to lie on a 
single curve. Thus, the dependence of bond stiffness on ordering, composition, or chemical system can be explained 
almost entirely in terms of bond length changes. Remarkably, other effects on force constant stiffness, such as changes 
in charge density associated with different configurations and chemical systems, are small. This may be because our 
study is limited to noble-metal intermetallics. 

For each pair type, we constructed a relationship between force constant stiffness and bond length. Our goal is 
to parameterize the stiffness of the stretching and bending force constants as a function of bond length, then use 
"fitted" force constants from these functions to predict vibrational entropy. Stretching terms were fit to a second- 
order polynomial; bending terms were fit to a linear function. The fit was restricted to first-nearest neighbors force 
constants. The rms fitting error over all stretching terms was 0.102 eV /A?\ the rms fitting error over all bending 
terms was 0.048 eV/A^. 

Although attempts to use longer-ranged fitted force constants were unsuccessful, the errors introduced by using only 
nearest-neighbor force constants was small. In the systems studied, vibrational entropy converged quickly with respect 
to neighbors: the difference between S^H) calculated using only first nearest-neighbors and Syib calculated using all 
neighbors was typically around ~0. 00-0. 02 fcs,with a maximum difference of 0.03 ks in LI2 PdsCu. Previous studies 
on intermetallics and group-IV semiconductors have also shown that vibrational entropy can converge quickly with 
respect to neighbors I2,i3,37_ 

The results of using fitted force constants to predict vibrational entropies are shown in Table 1. For all structures. 
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the vibrational entropies obtained using fitted force constants agree well with the vibrational entropies obtained 
using directly calculated force constants. The rms error between ^caic.oW neigh,, /c^ gfzt,inn,sb structures 
investigated was 0.032 fcs^with a maximum error of 0.058 ks in LI2 PdsAu. The errors are ~l-2% of the absolute 
vibrational entropy. These calculations show the predictive power of using fitted stretching-bending force constants 
to predict vibrational entropies. 

We also tested the accuracy of using fitted force constants, in representing the effect of homogeneous volume changes. 
In metals, pairwise expansion can usually not capture the effect of such electron density changes of the energy. The 
force constants and vibrational entropy of LI2 and SQS-8 structures for PdsCu, Pd^Au, CusPd, CusAu, and AusPd 
were recalculated at a volume ^2% larger than the equilibrium volmne. For the LI2 PdsCu structure, the vibrational 
entropies were also calculated at a 4% and 6% volume increase were also calculated. These direct results were then 
compared with predicted results obtained from using force constants derived from equilibrium volume calculations. 
None of the force constants at these increased volumes were included in the fit. The results of using these fitted force 
constants to predict vibrational entropy are shown in Table 2. For all structures at 2% increased volume, the rms 
error between and S^^*-!""'*'' was 0.033 ke, with a maximum error of 0.053 ke in LI2 PdgAu. 

Thus, even when the force constant data is not included in the fit, the length dependent force constant function gives 
force constants that accurately predict vibrational entropies. 

Often, one is interested in vibrational entropy differences, rather than absolute vibrational entropies. We used 
vibrational entropy data from Tables 1-2 to obtain two types of vibrational entropy differences that are of interest: 
the vibrational entropy difference between ordered and disordered phases, and vibrational entropy difference between 
a structure at two volumes separated by ~2%. The latter is related to the thermal expansion by the expression = 
where B is the bulk modulus. This data is shown in Table 3. In the systems studied, we estimated the effect 
of vibrations on calculated phase boundaries by using the following procedure: The change in transition temperature 

when vibrations are included in the phase diagram calculation is given by T^o2fig+vib ^ '^co^fig only [ ^ + as"^'' 

\ fon/,9 

Thus, ^ga^fi determines the effect of lattice vibrations on calculated phase diagrams. The configurational entropy 

config 

per atom for a binary solid solution depends on the state of short-range order. It has a maximum value for an ideal 
(random) solution of 8"^^^^^ = k[chic+ {1 — c) ln(l — c)]. Hence, the maximum value of ^3"^^-^ in the systems 
studied is 0.562 fcs/atom: this value occurs when a fully ordered state {S^^^f^g = 0) at c = 0.25 transforms to a fully 

random solid solution. In most cases, the actual value of ^S^^j-^ will be smaller than this, owing to short-range 
order in the disordered state and also some disorder in the low temperature phase. Nevertheless, our approximate 
estimate shows that phase boundaries in the systems studied would change by ~5-13% if vibrations were included. 
The data in Table 3 can be used to examine the ability of transferable force constants to predict vibrational entropy 

differences. The rms error between calculated and fit differences AS'^^j,''^' - AS'f^*' was 0.030 

ks, with a maximum error of 0.063 ks in CusAu. For vibrational entropy differences between structures at volumes 

2% apart, the rnis error between calculated and fit differences was 0.016 ks with a maximum error of 0.030 ks in LI2 
AusPd. Thus, the transferable force constants used here are able to predict vibrational entropy differences reasonably 
well. 

In examining the data, one can make a few observations with respect to the accuracy of using transferable force 
constants. First, errors in absolute entropies from using transferable force constants (Tables 1-2) were typically ~0.00- 
0.05 fcsfor all structures. This was true for structures with both large and small absolute vibrational entropies. Second, 
in all cases, the transferable force constants were able to correctly predict the relative hierarchy of entropy differences. 
Thus, large entropy differences are predicted to be large; small entropy differences are predicted to be small. A small 
entropy difference was never predicted to be large, or vice versa. Third, when using transferable force constants, errors 
in entropy differences between the same structure at different volumes tended to cancel; errors in entropy differences 
between different structures (at the same composition) did not tend to cancel as much (Table 3). Thus, it is likely 
that there are small structure-dependent contributions to the force constant stiffness that are not captured by our 
length-dependence force constant model. Fourth, using transferable force constants gave small absolute errors on 
vibrational entropy differences (~0. 00-0. 07 fcs). For small vibrational entropy differences, this sometimes resulted 
in ^15-20% errors, whereas the percentage errors in larger vibrational entropy differences were smaller. Lastly, this 
method can greatly reduce the computational cost of calculating vibrational properties of intermetallics. 



IV. CONCLUSION 



The vibrational entropies for 16 structures in the Ag-Au, Au-Cu, Au-Pd, Cu-Pd, and Cu-Pt systems have been 
calculated. A simplified model was used with only stretching and bending terms, making the force constant matrices 
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independent of symmetry. With these approximations, the form of the force constant matrices was independent of 
crystal symmetry. We found that the variation of force constants with ordering, composition, or chemical system can 
be explained almost entirely through changes in bond length. This method represents a promising way to include 
vibrational effects in phase diagram calculations at a moderate computational cost. 
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FIG. 1: Stiffness of bonding and stretching force constants vs. bond length for Au-Au bonds. Force constants for the same 
system are represented by the same symbol (for example, all Au-Au bonds for LI2 AusCu use the same symbol). Lines are 
drawn as guides to the eye. 
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FIG. 2: Stiffness of bending and stretching force constants vs. bond length for Cu-Cu bonds. Lines are drawn as guides to the 
eye. 
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FIG. 3: Stiffness of bending and stretching force constants vs. bond length for Pd-Pd bonds. Lines are drawn as guides to the 
eye. 




o 

§ 2.0 2.5 3.0 3.5 4.0 4.5 5.0 
Bond length (Angstroms) 



FIG. 4: Stiffness of bending and stretching force constants vs. bond length for Cu-Pd bonds. Lines are drawn as guides to the 
eye. 
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FIG. 5: Stiffness of bending and stretching force constants vs. bond length for Au-Pd bonds. Lines are drawn as guides to the 
eye. 
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FIG. 6: Stiffness of bending and stretching force constants vs. bond length for Au-Cu bonds. Lines are drawn as guides to the 
eye. 
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TABLE L Vibrational entropy data. Listed are vibrational entropies calculated with all neighbors and full force constant 
matrices, vibrational entropies calculated with all neighbors and stretching-bending force constants, calculated vibrational 
entropy with first-nearest neighbors and full force constant matrices, calculated vibrational entropy with first-neighbors and 

strctcliing-bcnding force constant matrices, and fit-vibrational entropies using first-neighbors and stretching-bending force 
constants. Also listed is(S'^"^'''"" — S^l^'^""), the error introduced by using fit force constants, the stretching-bending force 
constant model, and first-nearest neighbor force constants. Errors are 0.015 ks unless otherwise indicated. All numbers are 
rounded to 0.01 ks- 
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Structure 


Neighbors used in 
Ml force constant 
Calculation 


All neighbors, fiill 
force constants (^b) 


ri,lnn 

Inn, liill force 
constants 
«:b) 


Fit stretching and 
bending terms, Inn 
(« 




Pd3CuLl2(2%) 


5 


-4.41 


-4.44 


-4.42 


0.01 


PcbCuLl2(4.5%) 


5 


-4.30 


-4.33 


-4.31 


0.01 


PdjCuLl2(6.0%) 


5 


-4.14 


-4.17 


-4.16 


0.02 


PdjCu SQS8 (2%) 


2 


-4.39 


-4.37 


-4.40 


0.01 


PdbAuLl2(2%) 


3 


-4.08 


-4.06 


-4.02 


-0.05 


PdbAu SQS8 (2%) 


2 


-4.01 


-4.00 


-3.99 


-0.02 


Cu3PdLl2(2%) 


3 


-4.61 


-4.58 


-4.64 


4fl.03 


CujPd SQS8 (2%) 


2 


-4.63 


-4.62 


-4.62 


-0.01 


CujAu LI2 (2%) 




-4 18 


-4 15 


-4 13 


-0 05 


CusAu SQS8 (2%) 


2 


-4.12 


-4.14 


-4.13 


-tO.Ol 


Au3PdLl2(2%) 


5 


-3.24 


-3.26 


-3.27 


-tO.03 


AuaPd SQS8 (2%) 


2 


-3.28 


-3.27 


-3.23 


-0.05 



'vib > 



TABLE II: Vibrational entropies at increased volume. Percentage increase in volume is indicated in parenthesis. Listed are 

vibrational entropies calculated with all neighbors and full force constant matrices, vibrational entropies calculated with fit 
first-nearest neighbor stretching-bending force constants, and (5'^^^'='"" _ sf^ll'^""'"''). All numbers are rounded to 
0.01 kB. 
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System 


A c t»/c , aU neighbors , JUll fc 


A C Jit,lnn,sb 


A ^ cede , all neigh, full fc _ fit,\tin,sb 
vib vib 




AU neighbors, full force 


Fit stretching and 








constants 


bending 








(fe) 


terms, Inn(^B) 






PdsCuSQSS-PdjCuLh 


+0.03 


+0.05 




+0.02 


PdsAuSQSS-PdsAuLli 


+0.07 


+0.02 




+0.05 


Cu3PdSQS8-Cu3PdLl2 


+0.03 


+0.00 




+0.03 


CujAu SQS8 - CujAu Lh 


+0.07 


+0.00 




-0O6 


Au3PdSQS8-Au3PdL12 


-0.05 


-0.04 




+0.01 


Pd3Cu LI2 (2%) - Pd3Cu LI2 


+0.14 


+0.15 




+0.00 


PdjCu SQS8 (2%) - PdjCu SQS8 


+0.14 


+0.15 




-001 


PdsAu LI2 (2%) - PdsAu Lh 


+0.15 


+0.15 




+0.01 


Pd3Au SQS8 (2%) - Pd3Au SQS8 


+0.15 


+0.16 




-002 


Cu3PdLl2(2%)- CujPdLh 


+0.12 


+0.12 




+0.00 


CusPd SQS8 (2%) - CusPd SQS8 


+0.13 


+0.14 




-001 


Cu3AuLl2(2%)- CU3AULI2 


+0.12 


+0.15 




-002 


CU3AU SQS8 (2%) - CiijAu SQS8 


+0.12 


+0.14 




-0.02 


AusPd LI2 (2%) - AujPd LI2 


+0.16 


+0.13 




-0O3 


AusPd SQS8 (2%) - AujPd SQS8 


+0.17 


+0.17 




+0.00 


Pd3CuLl2(4.5%) - PdjCuLh 


+0.25 


+0.26 




+0.00 


PdjCu LI2 (6.0%) - PdsCu LI2 


+0.42 


+0.41 




+0.00 



TABLE III: Vibrational entropy differences. Listed are vibrational entropy differences upon disordering and vibrational entropy 
differences between the same structure at different volumes. Percentage increase in volume is indicated in parenthesis. Also 
shown are errors in differences, AS^fJ'^'"""'^*^'''""-^ - A5^^^*'i"". 



